function normMod(cond_one,cond_two,density,nSurrogates,outFile)


%density=0.05:0.05:0.40;

%[th_min,th_one,th_two]=findRightThoSubjs(cond_one,cond_two,th)

[nEdge_one,connected_one,obsCixD_one,obsQxD_one,avgSurrQxD_one,stdSurrQxD_one,obsQsdxD_one]=normalize(cond_one,density,nSurrogates);
[nEdge_two,connected_two,obsCixD_two,obsQxD_two,avgSurrQxD_two,stdSurrQxD_two,obsQsdxD_two]=normalize(cond_two,density,nSurrogates);

save(outFile,'density','nSurrogates','nEdge_one','connected_one','obsCixD_one','obsQxD_one','avgSurrQxD_one','stdSurrQxD_one',...
'nEdge_two','connected_two','obsCixD_two','obsQxD_two','avgSurrQxD_two','stdSurrQxD_two','obsQsdxD_one','obsQsdxD_two');




function [nEdge,connected,obsCixD,obsQxD,avgSurrQxD,stdSurrQxD,obsQsdxD]=normalize(cond,density,nSurrogates)

nDensity=length(density);
nSubjs=size(cond,1);
nch=size(cond,2);
obsCixD=zeros(nDensity,nSubjs,nch);
obsQxD=zeros(nDensity,nSubjs);
obsQsdxD=zeros(nDensity,nSubjs);
avgSurrQxD=zeros(nDensity,nSubjs);
stdSurrQxD=zeros(nDensity,nSubjs);

connected=zeros(nDensity,nSubjs);
nEdge=zeros(nDensity,nSubjs);
for j=1:nDensity
    for i=1:nSubjs
        
        [thM,connected(j,i),nEdge(j,i)]=isConnected(squeeze(cond(i,:,:)),density(j));    
        [obsCixD(j,i,:),obsQxD(j,i),avgSurrQxD(j,i),stdSurrQxD(j,i),obsQsdxD(j,i)]=surrogate(thM,nSurrogates);      
    end
end


function [obsCi,obsQ,avgSurrQ,stdSurrQ,obsQsd]=surrogate(M,nSurrogates)

%nSurrogates=50;
rewiredP=1;


nch=length(M);
[obsCi,obsQ,obsQsd]=computeModMeasures(M);
surrogatesQ=zeros(nSurrogates,1);
surrogatesCi=zeros(nSurrogates,nch);
for i=1:nSurrogates
        rM=randmio_und_connected(M,rewiredP);
        %rM=randomizer_bin_und(M,rewiredP);
        [surrogatesCi(i,:),surrogatesQ(i)]=computeModMeasuresSURR(rM);
end
avgSurrQ=mean(surrogatesQ);
stdSurrQ=std(surrogatesQ,[],1);




function [thM,connected,nEdge]=isConnected(M,density)
nch=length(M); 
thM=threshold_proportional(M,density);
thM(thM>0)=1;
thM=tril(thM,-1);
thM=thM+thM';
if(isequal(thM,thM.')) 
    [disc,pred,cl] =graphtraverse(sparse(thM),1,'METHOD','DFS','DIRECTED',false);
    connected=0;
    if(length(disc)==nch)
            connected=1;
    end
    nEdge=sum(sum(tril(thM,-1)));  
else
    error('not symmetric matrix after threshold!!!')
end

function [Ci,Q,Qsd]=computeModMeasures(M)
ntimes=100;
s=size(M);
Ci=zeros(1,s(1));
%PartCoef=zeros(1,s(1));
%Zintra=zeros(1,s(1));
Q=zeros(1,ntimes);
for i=1:ntimes
[Ci,Q(i)]=modularity_louvain_und(M);

%PartCoef=participation_coef(M,Ci)';
%Zintra=module_degree_zscore(M,Ci)';
end
Qsd=std(Q);
Q=mean(Q);

function [Ci,Q]=computeModMeasuresSURR(M)

s=size(M);
Ci=zeros(1,s(1));
%PartCoef=zeros(1,s(1));
%Zintra=zeros(1,s(1));

[Ci,Q]=modularity_louvain_und(M);

%PartCoef=participation_coef(M,Ci)';
%Zintra=module_degree_zscore(M,Ci)';

